# Read ImageJ Result files that contain features of each WPB
features <- NULL
features <- read.table("/path/*_Results_Features.csv", header=T, sep=",") # e.g. A2_Results_Features.csv

# Read files that match a WPB to a certain cell
WPBcellID <- NULL
WPBcellID <- read.table("/path/*_Results_WBP_cell_ID.csv", header=T, sep=",") # e.g. A2_Results_WBP_cell_ID.csv

names(WPBcellID) <- c("no","cell.ID","Slice")
cellID <- WPBcellID$cell.ID

# Attach the column to the Features that matches each WPB with a cell
features <- cbind(features,cellID)

attach(features)

# Add plateID to the table as a column
plateClearBack <- sub(":.*", "", Label)
plateID <- sub("^.{13}", "", plateClearBack)
features <- cbind(features,plateID)	

# Add row number to the table as a column
wellrowClearFront <- sub("^.{18}", "", Label) # Delete the string before the row (A-H) identifier 1-8
wellrowID <- as.numeric(sub("0.*", "", wellrowClearFront)) # Delete the string with the 0 after the row (A-H) identifier 1-8
features <- cbind(features,wellrowID)

# Add column number to the table as a column
wellcolClearFront <- sub("^.{20}", "", Label) # Delete the string before the column identifier 1-12
wellcolID <- as.numeric(sub("_.*", "", wellcolClearFront)) # Delete the string with the _ after the column identifier 1-12 000 and convert to numeric from character
wellcolID <- wellcolID / 1000	# Remove the 000
features <- cbind(features,wellcolID)

# Add the field of view number to the table as a column
fieldID <- sub("^.{32}", "", Label) # Delete the string before the field identifier 1-5
features <- cbind(features,fieldID)

attach(features)

# One value per well calculation for one plate. 

# (1) 1_PercentWPBarea: area ratio of WPBs with Feret’s diameter larger than 1.5um 
# (2) 2_TotalWPBnrPerCellNR: WPB number per cell
# (3) 3_TotalWPBrawIntDenPerCellNR: WPB intensity per cell
# (4) 4_TotalFoVpercentWPBnr:	ratio of WPB number with Feret’s diameter larger than 1.5um

PercentWPBarea <- NULL # for (1)
TotalWPBnrPerCellNR <- NULL # for (2)
TotalWPBrawIntDenPerCellNR <- NULL # for (3)
TotalFoVpercentWBPnr <- NULL # for (4)

n=0
for (r in 1:8) { 
 for (c in 1:12) {

   n=n+1

   FoVpercentWBParea <- NULL  # for (1)
   FoVTotalWPBnrPerCellNR <- NULL  # for (2)
   FoVTotalWPBrawIntDenPerCellNR <- NULL  # for (3)
   FoVpercentWBPnr <- NULL # for (4)

   for (f in 1:5) {
   FoVfeature <- features[wellrowID==r & wellcolID==c & fieldID==f & cellID>0,] 
   FoVcellNR <- length(levels(as.factor(FoVfeature$cellID))) # FoVfeature$cellID Number of cells per image (field of view (FoV)) 
   FoVTotalWPBareaPerCellNR <- sum(FoVfeature$Area)/FoVcellNR	# Total area of WPBs per cell for (1)

   over15umfeature <- features[wellrowID==r & wellcolID==c & fieldID==f & cellID>0 & Feret>1.5,] 
   FoVover15umWPBarea <- sum(over15umfeature$Area)/FoVcellNR	# Total area of WPBs longer than 1.5um Feret per cell for (1)
   FoVpercentWBParea[f] <- FoVover15umWPBarea/FoVTotalWPBareaPerCellNR # for (1)
    if(!is.finite(FoVpercentWBParea[f])||FoVpercentWBParea[f]==1) FoVpercentWBParea[f] = 0 # avoid DiV/0 (=no cell) or 1 WPB

   FoVTotalWPBnr <- length(FoVfeature$RawIntDen)   # Total NR of WPB in current FoV for (2)
   FoVTotalWPBnrPerCellNR[f] <- FoVTotalWPBnr/FoVcellNR # for (2)
    if(!is.finite(FoVTotalWPBnrPerCellNR[f])||FoVTotalWPBnrPerCellNR[f]==1) FoVTotalWPBnrPerCellNR[f] = 0 # avoid DiV/0 (=no cell) or 1 WPB
    
    FoVpercentWBPnr[f] <- length(over15umfeature$Area)/FoVTotalWPBnr # (4) Ratio: number of WPBs longer than 1.5um Feret per total nr of WPBs (per FoV)
    if(!is.finite(FoVpercentWBPnr[f])||FoVpercentWBPnr[f]==1) FoVpercentWBPnr[f] = 0 # avoid DiV/0 (=no cell) or 1 WPB

   FoVTotalWPBrawIntDen <- sum(FoVfeature$RawIntDen)   # Sum of Total Intensities of WPB in current FoV for (3)
   FoVTotalWPBrawIntDenPerCellNR[f] <- FoVTotalWPBrawIntDen/(FoVcellNR * 1000000) # for (3) 
    if(!is.finite(FoVTotalWPBrawIntDenPerCellNR[f])) FoVTotalWPBrawIntDenPerCellNR[f] = 0 # avoid DiV/0
    
    
   }

  PercentWPBarea[n] <- mean(FoVpercentWBParea) # (1)
  TotalWPBnrPerCellNR[n] <- mean(FoVTotalWPBnrPerCellNR) # (2)
  TotalWPBrawIntDenPerCellNR[n] <- mean(FoVTotalWPBrawIntDenPerCellNR)  # (3)
  TotalFoVpercentWBPnr[n] <- mean(FoVpercentWBPnr) # (4)
 }
}

output_wellID <- c("A01","A02","A03","A04","A05","A06","A07","A08","A09","A10","A11","A12","B01","B02","B03","B04","B05","B06","B07","B08","B09","B10","B11","B12","C01","C02","C03","C04","C05","C06","C07","C08","C09","C10","C11","C12","D01","D02","D03","D04","D05","D06","D07","D08","D09","D10","D11","D12","E01","E02","E03","E04","E05","E06","E07","E08","E09","E10","E11","E12","F01","F02","F03","F04","F05","F06","F07","F08","F09","F10","F11","F12","G01","G02","G03","G04","G05","G06","G07","G08","G09","G10","G11","G12","H01","H02","H03","H04","H05","H06","H07","H08","H09","H10","H11","H12")

#-- Save (1) PercentWPBarea
output <- NULL
for (a in 1:96) {output[a] <- c("plateID")}	# Add the column with plate IDs

output <- cbind(output,output_wellID)	# Add the column with the well IDs

output <- cbind(output,PercentWPBarea)	# Add the column with the values

output <- as.data.frame(output)

write.table(output, file = "/path/1_PercentWPBarea/plateID.txt", sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE) 

#-- Save (2) TotalWPBnrPerCellNR
output <- NULL
for (a in 1:96) {output[a] <- c("plateID")}	# Add the column with plate IDs

output <- cbind(output,output_wellID)	# Add the column with the well IDs

output <- cbind(output,TotalWPBnrPerCellNR)	# Add the column with the values

output <- as.data.frame(output)

write.table(output, file = "/path/2_TotalWPBnrPerCellNR/plateID.txt", sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE) 

#-- Save (3) TotalWPBrawIntDenPerCellNR
output <- NULL
for (a in 1:96) {output[a] <- c("plateID")}	# Add the column with plate IDs

output <- cbind(output,output_wellID)	# Add the column with the well IDs

output <- cbind(output,TotalWPBrawIntDenPerCellNR)	# Add the column with the values

output <- as.data.frame(output)

write.table(output, file = "/path/3_TotalWPBrawIntDenPerCellNR/plateID.txt", sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE) 

#-- Save (4) TotalFoVpercentWBPnr
output <- NULL
for (a in 1:96) {output[a] <- c("plateID")}	# Add the column with plate IDs

output <- cbind(output,output_wellID)	# Add the column with the well IDs

output <- cbind(output,TotalFoVpercentWBPnr)	# Add the column with the values

output <- as.data.frame(output)

write.table(output, file = "/path/4_TotalFoVpercentWBPnr/plateID.txt", sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE) 
